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Landslides are a common type of disaster in Indonesia, especially in steep- 
slope areas. The landslide process can be well understood by measuring the 
surface deformation. Currently, there are no practical solutions for measuring 
surface deformation at landslide locations other than field surveys in the 
Pacitan Regency. We apply LiCSBAS, to identify surface deformation in 
several landslide locations in a specific non-urban area with mixed 
topographical features. LiCSBAS is a module that utilizes data from the 
project of looking inside the continent from space (LiCS), using the new small 
baseline area subset (NSBAS) method. This study utilizes the leaf area index 
(LAD to validate the ability of LiCSBAS to detect surface deformation values 
at landslide locations. The study succeeded in identifying surface 
deformations at 100 landslide locations, with deformation values ranging from 
15.1 to 10.9 millimeters per year. Most of the landslide locations are closely 
related to volcanic rocks and volcanic sediments on slopes of 30-35°. The 


NSBAS method in the LiCSBAS module can reduce gaps error in the 
sentinel-1 image network. However, the utilization of the C-band at a pixel 
size of 100 meters made surface deformation only well detectable in a large 
open landslide area. 
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1. INTRODUCTION 

Soil surface deformation is a dynamic process on the earth's surface [1]. This process occurs naturally 
or due to human intervention [2]. One form of such deformation is a landslide. Landslide is a mass movement 
of soil or rock on the disrupted slope [3]. Landslides have caused casualties and material losses in many 
countries [4]. Pacitan Regency is the area with the highest number of landslide occurrences in East Java 
Province. The worst landslide events occurred in 2017, along with the Cempaka tropical cyclone; more than 
210 landslides were causing 19 deaths, destroying 615 houses, and causing loss of more than 615 billion rupiahs 
[5]. Therefore, strategic planning to reduce losses and casualties is necessarily conducted through data 
inventory (mapping) and ground movements identification [6]—[9]. 

Methods for detecting and monitoring ground movements are continuously evolving using remote 
sensing imagery, both optical and synthetics aperture radar (SAR) and field measurement data [2]. Some 
advantages of using satellite monitoring include having a large coverage area, immediate data acquisition, and 
ultimately reducing operational costs. One of them is using differential interferometry synthetic aperture radar 
(DInSAR). DInSAR works by utilizing different phases of SAR images in the same area at different times [10]. 
The use of DInSAR for the detection and monitoring of ground motion includes the persistent scatterer 


Journal homepage: http://ijece.iaescore.com 


Int J Elec & Comp Eng ISSN: 2088-8708 O 2615 


interferometry synthetics aperture radar (PSInSAR) and small baseline area subset (SBAS) method [2], 
[11]-[15]. PSInSAR works using a single master interferogram to form persistent scatterer (PS) [16]. 
Meanwhile, SBAS works by utilizing multiple master interferograms with smaller temporal and spatial 
baselines [17]. The study of surface deformation in non-urban areas usually uses the SBAS method by utilizing 
several master interferograms [12], [14], [18]. Moreover, SBAS does not depend on PS density to evaluate 
surface deformation [18]. However, the SBAS method uses large amounts of data and requires compatible 
devices for processing on a long temporal scale and regional-level analysis. To overcome these limitations, 
users can use LiCSBAS, an open-source InSAR multi-temporal analysis module [19]—[21]. LiCSBAS utilizes 
sentinel-1 SAR data that LiCSAR has processed, so users do not need to download enormous datasets on their 
computers [20]. LiCSBAS uses the Python3 and the bourne again shell (Bash) to write commands. LiCSBAS 
allowed us to measure surface deformation at a millimeter-scale [19]. 

The current use of LiCSBAS focuses on large-scale surface deformation in tectonic and volcanic 
studies [21]. LiCSBAS, on the other hand, has the potential to be used in hydrosphere, cryosphere, and mass 
movement studies [20]. Unfortunately, there are no practical solutions for addressing surface deformation in 
landslide locations with various physical and environmental conditions other than field surveys. Another 
challenging issue is related to the implementation of LiCSBAS in a non-urban area with different lithological 
and topographical characteristics. Therefore, this study aims to implement the LiICSBAS module on a multi- 
temporal InSAR to detect surface deformation values at landslides locations in the Pacitan Regency, which has 
the characteristics described above. This study also utilizes overlay analysis to determine the distribution of 
landslide locations based on lithological formations and slopes. In addition, this study uses the Leaf Area Index 
(LAD) to determine the vegetation dynamics that can affect the coherence of the interferogram. 


2. RESEARCH METHOD 
2.1. Study area 

The study area of Pacitan Regency, East Java Province, Indonesia, consists of 12 sub-districts with 
about 1,417 km? and occupied by 586,110 people in 2020 [22]. Land-use/land-cover data from the Indonesian 
geospatial information agency (BIG) shows that the Pacitan Regency is dominated by shrubs, covering more 
than 36,000 hectares. In contrast, others types are dryland farming (31,444 hectares), mixed plantations 
(18,540 hectares), rainfed rice fields (18,540 hectares), and settlement areas (15,031 hectares) [23]. This study 
used 358 landslide locations. As many as 258 landslide locations are the data from the regional disaster 
management agency (BPBD) of Pacitan Regency inventory in 2017-2020. The remaining 100 landslide 
locations are the results of identification through field surveys and interpretation of high-resolution Airbus and 
the Centre National d'Etudes Spatiales (CNES) imagery in 2017-2020 accessed via Google Earth. 

Hilly topography dominates the study area, which is also part of the Southern Mountains [24]. The 
downstream of the Grindulu River and the Lorok River have a relatively flat topography. Table 1 shows the 
results of slope classification in Pacitan Regency. The central part of the Pacitan Regency has a steeper slope 
compared to the southern part. However, as indicated in Figure 1, numerous high-degree slopes were also 
found in the south-eastern part of the area, particularly in the northern part of the Ngadirojo Sub-District. The 
Arjosari Formation is the widest, and the Kalipucang Formation is the narrowest, according to the area of each 
lithological formation as shown in Figure 2 and Table 2. 


Table 1. Percentage areas of slope classes 
No Slope Area (hectares) (percent) 


1 0-5° 13,250 9.3 
2 5-10° 25,755 18.1 
3 10-15° 30,399 21.4 
4 15-20° 27,744 19.5 
5 20-25° 21,056 14.8 
6 25-30° 13,573 9.5 
7 >30° 9,951 7 


Source: Calculated from Digital Elevation 
Model-National (DEMNAS) [25] based on [26] 


On the other hand, leaf area index (LAI) indicates the density of plant canopy in a particular area and 
is usually defined in m? [27], [28]. Figure 3 illustrates the average value of LAI in the Pacitan Regency from 
2017-2020 predicted from sentinel-2A Multi Spectral Instrument (MSI) imagery processed using the Google 
Earth Engine (GEE). The lowest LAI values are illustrated as red-colored zones along the Grindulu River and 
the Lorok Rivers, indicating sedimentary material and sparse vegetation cover. Meanwhile, high LAI values 
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are illustrated as blue-colored zones in the central and eastern parts of the Pacitan Regency with hilly 
topography and covered with agricultural and non-agricultural vegetation. 
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Figure 1. Slope class distribution in Pacitan Regency, calculated from DEMNAS [25] based on [26] 
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Figure 2. Lithological formations in Pacitan Regency [29], [30] 


Table 2. Distribution of rock formation 


No Formation Rock materials Area Area 
(hectares) (%) 
1 Alluvial alluvium 7,397 5.2 
2 Arjosari polymictic breccias, sandstone, and conglomerates alternating with volcanic 33,150 23.3 
rock 
3 Dacite dacite, alluvium, and another intrusive rock 1,101 0.7 
4 Dayakan sandstone and claystone 1,543 1 
5 Intrusive Rock intrusive rock 3,250 22 
6 Jaten conglomerate, sandstone, mudstone, lignite, shale, and tuff 7,743 5.4 
7 Kalipucang conglomerates and clay 18 0.01 
8 Mandalika volcanic breccias and lava tuffs altering with sandstone and siltstone 26,273 18.5 
9 Nampol sandstone, siltstone, limestone, and lignite with conglomerate and breccia 3,624 2.5 
10 Oyo sandstone, siltstone, limestone, and marl 4,466 3.1 
11 Semilir tuff, sandstone breccia, and siltstone 2,974 2 
12 Watupatok basalt pillow lava, sandstone, claystone, and chert 14,676 10.3 
13 Wonosari reef limestone, limestone deposits, sandy limestone, and marl 30,114 21.2 
14 Wuni volcanic breccias, tuffs, and sandstone alternating with lignite and limestone 5,380 3:7 


Source: Lithological formations of Pacitan Regency [29], [30] 
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Figure 3. Average LAI from 2017 to 2020, calculated from sentinel-2A multi spectral instrament (MSI) 
imagery [31] 


2.2. Method of SAR data processing 

The LiCSBAS module utilizes the LiCSAR product, so users do not need to create interferograms 
from SLC data on the local computer [19]. This study utilizes 190 Sentinel-1 images of ascending 
Interferometric Wide, single look complex, and vertical-vertical polarisation (IW SLC VV). The images are 
automatically processed in LiCSAR to form 812 interferograms with a resolution of 0.0001 or 100 meters [19], 
[20]. The 190 sentinel-1 images have a 3.24-year observation period, starting from 07-05-2017 and terminating 
on 01-08-2020, corresponding to the period included in the landslide data inventory. LiCSBAS, coded in 
python and bourne again shell (bash), was conducted in this study using Ubuntu 20.04 long term support (LTS) 
software [19]. In general, there are two categories of data processing: data preparation and multi-temporal 
analysis [19]. In the research [19], [20] for detailed data processing procedures. First, download unwrapped 
phase and coherence GeoTIFF data from the Centre for the Observation and Modelling of Earthquakes, 
Volcanoes and Tectonics (COMET)-LiCS website, www.comet.nerc.ac.uk [19]. Second, convert geographic 
tagged image file format (GeoTIFF) files to a single-precision floating-point format, then used the generic 
atmospheric correction online service (GACOS) for InSAR data to correct tropospheric noise [19]. 
Furthermore, mask the de-correlated parts of the interferogram where the phase noise coherence estimates are 
0.5 [19], and clip the rectangular area of interest at-7.89° (North), -8.28° (South), 110.90° (West), 111.46° 
(East). 

There are four stages of multi-temporal analysis. First, identify the interferograms based on coherence 
and unwrapped data coverage that may experience spatial and temporal de-correlation. The next step is by 
checking the loop closure phase to remove interferograms with many unwrapping errors [19]. The (1) describes 
the loop closure [32]. Ø1, Ø2, Ø; are images, and Ø12, G23, Ø13 are the unwrapped interferograms. The loop 
phase usually has a value close to zero if there is no unwrapping error in the three interferograms [19]; 


i23 = O12 + O23 — Dis (1) 


Second, to obtain the time-series velocity displacement as line of sight (LOS) displacement, invert the stack of 
unwrapped data [19]. According to [19] and [33], the time series velocity calculation utilizes the NSBAS 
method by inverting the interferogram network, as shown in (3). The NSBAS technique is a modification of 
the Small Baseline approach in (2) to overcome the limitation of a time gap on the network [19], [34]. 
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Third, calculate the standard deviation of the velocity, mask the noisy pixels, and filter the time series 
velocity [19]. Fourth, calculate and visualize the value of surface deformation at each landslide location. The 
surface deformation value or LOS displacement value obtained from the interferogram is the result of 
measuring the movement between the radar antenna and a location on the earth's surface [15]. This value is 
represented by a size of 100 by 100 meter pixels resolution. 


3. RESULTS AND DISCUSSION 
3.1. Surface deformation at landslide locations 

Surface displacement in the study area varies from -34 to 25 millimeters per year in an area of 33,831 
hectares, or 23.8 percent of the Pacitan Regency. As illustrated in Figure 4, the type of deformation observed 
was subsidence and uplift, covering 10,062 hectares (7.1%) and 23,768 hectares (16.7%), respectively. The 
most extensive deformation was found in the agricultural and bare-land areas, totaling 21,272 hectares, 
followed by a built-up area (7,457 hectares) and non-cultivated land (4,367 hectares). 

Of the 358 landslide locations, 46 landslide locations have negative displacement values ranging from 
15.1 millimeters per year to 0.05 millimeters per year. On the other hand, 54 landslide locations had positive 
displacement values ranging from 0.3 to 10.9 millimeters per year. Meanwhile, the displacement value at 258 
landslide locations remains undetected. This undetected displacement at landslide locations is associated with 
a small landslide area covered with vegetation, as shown in Figure 5. 

Surface deformation is undetected in most areas of the Pacitan Regency due to the low coherence 
values found in several areas. This circumstance might be related to the geometric distortion effect due to the 
side-looking direction of the image acquisition in the SAR system [35]. For example, the shadow can occur 
when hills block the radar beam so that the sensor receives no backscatter from the object, which is 
characterized by dark areas in the coherence image [36]. Furthermore, dense vegetation might cause the 
backscatter to be weaker, reducing the coherence value in the interferogram [21]. Therefore, the LiCSBAS 
module will automatically mask pixels with a coherence value of less than 0.5 [19]. 

High hilly areas covered with hardwood trees and dense canopy in the central part of the Pacitan 
Regency have low density of surface deformation pixel. The dominance of high LAI values indicates this 
circumstance. Figure 4 shows the distribution and pixel values of the surface deformation. The built-up region 
(in the red box) has a high surface deformation pixel density because the C-band SAR signal produces a high 
coherence value. The red boxes on the west and east sides represent the downstream of the Grindulu River in 
the Pacitan Sub-District and the downstream of the Lorok River in Ngadirojo Sub-District, respectively. In the 
west part of the Pacitan Regency, in the Sub-Districts of Donorojo, Punung, and Pringkuku, medium-density 
deformation pixels may be detected. Most of this area is open land with sparse vegetation. 

The spatial information of landslide zones is essential to analyze the LOS displacement value. For this 
purpose, landslides can be divided into depletion and accumulation zone [2], [37]. The depletion zone is a 
source of landslides characterized by a scarp at the top of the landslide [2], [26]. The accumulation zone is 
considered the fill zone of landslide material from the depletion zone [37]. The depletion zone has a negative 
LOS displacement value because the distance between a point on the earth's surface and the satellite radar 
sensor gets further away and vice versa in the accumulation zone [38]. 
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Figure 4. Distributions of surface deformation velocity pixels Figure 5. Small landslides on the side of 
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Figures 6 illustrates three landslide locations whose deformation from LOS displacement values were 
detected, both in the depletion and accumulation zones or one of them. The delineation lines on the scarp, 
depletion, and accumulation zones result from field surveys and identification of satellite images. The three 
landslide locations are included in the Arjosari and Mandalika Formation, with a slope class of 10-20°, each 
area of more than 1 hectare. Figure 6(a) is a landslide in Kalikuning Village, Tulakan Sub-District, with an 
elevation at point of interest (POI) 1 of 336 meters. POI 1 is located at the bottom of the landslide or 
accumulation zone, evidenced by a positive LOS displacement value of 9.4 millimeters per year. Figure 6(b) 
is a landslide in Kedungbendo Village, Arjosari Sub-District, with an elevation of POI 2 of 57 meters and POI 
3 of 245 meters. POI 2 is located in the accumulation zone, evidenced by a positive LOS displacement value 
of 8.8 millimeters per year. POI 3 is located in the depletion zone, with a negative LOS displacement value of 
-5.1 millimeters per year. Figure 6(c) is a landslide in Ponggok Village, Pacitan Sub-District, with an elevation 
at POI 4 of 345 meters. POI 4 is located in the depletion zone, evidenced by a negative LOS displacement 
value of -9.8 millimeters per year. 

In Figure 7, the LOS displacement time-series at each POI shows a trend indicating surface 
deformation still happens at landslide locations during the study period. In Figures 7(a) and 7(b), the trend line 
at POI 1 and POI 2 has an increasing trend with a positive average LOS displacement. The positive LOS 
displacement value indicates that the distance between the earth's surface and the satellite sensor is getting 
closer due to the increasing pile of material from the depletion zone. The opposite occurs in Figures 7(c) and 
7(d), where the average LOS displacement values at POI 3 and POI 4 are decreasing, which indicates the 
distance between the earth's surface and the satellite sensor is getting farther away caused by the decrease in 
slope materials towards the accumulation zone. 
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Figure 6. Detected LOS displacement at landslide zone and point of interest locations (a) detected 
accumulation zone in Kalikuning Village, Tulakan Sub-District, (b) detected depletion and accumulation 
zone in Kedungbendo Village, Arjosari Sub-District, and (c) detected depletion and accumulation zone in 

Ponggok Village, Pacitan Sub-District 
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In order to classify the activeness of landslides based on their surface deformation values, this study 
uses a threshold with a value of LOS displacement of 1.5 millimeters per year. According to [12], landslides 
are considered moving or active when at least one displacement pixel is in the landslide location with a value 
of LOS displacement of less than -1.5 or more than 1.5 millimeters per year. The NSBAS technique in the 
LiCSBAS module can identify the LOS displacement values at 100 landslide locations in Pacitan Regency. Of 
the 74 landslide locations that are still active, it consists of 30 landslide locations with negative displacement 
values between-1.5 and 15.1 millimeters per year and 44 other locations with positive values of 1.5 to 10.9 
millimeters per year. The other 258 locations could not be classified because the displacement values were not 
detected. According to [12], this is due to the low coherence value in most areas because of vegetation cover 
and shadow errors in hilly areas. 
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Figure 7. LOS displacement time-series for selected point of interest at the accumulation zone (a) POI 1, 
(b) POI 2; and at the depletion zone, (c) POI 3, and (d) POI 4 


3.2. Distribution of landslide locations based on lithological formations and slopes 

The overlay analysis between the geological map sheets of Pacitan and Ponorogo with the landslide 
location shows that the Arjosari formation has the highest number of landslide occurrences (191 locations; 
45 percent), followed by Mandalika formation (53 locations; 18.5 percent) and Jaten formation (36 locations; 
10 percent) as shown in Figure 8(a). The three rock formations are volcanic rock and old volcanic sediments 
that form steep slopes and dominate, especially in the eastern part of the Pacitan Regency [40]. These 
sedimentary rocks have a quickly loose and weathered structure; therefore, they have higher landslide 
susceptibility [41]. 

The formations with the least number of landslide occurrences are Wonosari (17 locations; 
4.7 percent), Intrusive Rock (16 locations; 4.4 percent), Wuni (14 locations; 3.9 percent), Nampol (12 locations; 
3.3 percent), Watupatok (11 locations; 3 percent), Semilir and Oyo (7 locations; 1.9 percent), Dacite and 
Dayakan (1 point; 0.2 percent). Meanwhile, there was no landslide during 2017-2020 recorded in the 
Kalipucang formation. Only a few landslides appeared in the western part of the Pacitan Regency, where the 
Wonosari Formation dominates. This formation forms low labyrinth-cone karst due to limestone weathering 
[40], [42]. Although the Wonosari formation occupies 30,114 hectares or 21.2 percent overall, only 17 landslide 
locations or 4.7 percent are associated with it, which means that the Wonosari Formation has lower landslide 
susceptibility. 

Based on the overlay results, it is known that most of the landslide’s location at 15-30, as shown in 
Figure 8(b). In more detail, the most landslide locations were found in the slope class 15-20 (82 locations; 
23 percent), followed by the slope class 20-25 (76 locations; 21%), class 25-30 (61 locations; 17%), class 
10-15 (53 locations; 14%), class 5-10° (37 locations; 10%), class 30-35° (25 locations; 7%), class 25-40 
(12 locations; 3%), class 0-5 (10 locations; 3 %), and class 40-45 (2 locations; 0.5%). 
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Based on the overlay between the slopes and the landslide locations, most of the landslides were 
recorded in the 15-20 slope class of 83 locations or 22.9 percent, followed by the 20-25 slope class of 76 
locations or 21.2 percent, and the 25-30° slope class of 61 locations or 17 percent. Based on Figure 8(b), the 
frequency of landslide locations increases until the slope class of 15-20, then decreases as the slope class 
increases. According to [43]-[45], the decrease in the number of landslides after the 15-20 slope class has a 
relationship with the high rate of erosion along with the increase in the slope value. As a result, the soil material 
cannot thicken, and stable hard rock is exposed. Meanwhile, erosion on riverbanks causes landslides on slopes 
of 0-5. 

Apart from natural conditions, human activities are also influencing landslide occurrences [46]. For 
example, surface modification by cutting slopes for road and house construction can reduce soil shear 
resistance; therefore, it cannot maintain slope stability under saturated conditions [47], [48]. As a result, The 
shear strength on the slope will increase while the resisting strength will decrease, thereby increasing the 
landslide probability [48]. 

Land cover type might also affect the occurrence of landslides. In this study, 66 landslide locations, 
or 18.4 percent, were found in built-up land and 156 locations, or 43.5 percent, in cultivated and open land. 
Both types of land cover indicate the intervention of human activities. Meanwhile, only 136 landslide locations, 
or 37.9 percent, were found on non-cultivated land, including forest. According to [49], naturally vegetated 
areas such as forests can withstand shear stress soil due to strengthening by roots. Furthermore, it can reduce 
the chance of landslides. The study using the combination of NSBAS and superimposed lithological formation 
and slopes can identify the distribution of landslide locations. Most of the landslide locations are located on 
volcanic rocks and volcanic sediments with slopes of 15-30°. The frequency of landslide occurrences increases 
until the slope class of 15-20°, then decreases as the slope class increases. 
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Figure 8. Frequency distribution of landslide in each (a) lithological formation and (b) slope class 


3.3. The advantages and limitations of using LiCSBAS module 

Sentinel-1 SAR data processing using LiCSBAS has several advantages. First, LiCSBAS can be a 
solution for mapping and assessing surface deformation velocity over a large area [19]. Second, the critical 
Sentinel-1 SAR data is accessible and already processed by LiCSAR [19], [20]. Third, users can freely access 
LiCSAR products through the COMET-LiCS website without downloading large datasets to their computers 
[19], [20]. The utilization of LiCSAR data can reduce processing time and memory usage. In addition, the use 
of Sentinel-1 SAR data allows temporal analysis possible every six days in the study area. Utilization of the 
NSBAS method in the LiCSBAS module can reduce the gap and the de-correlation between time in the network 
[19], [34]. On the other hand, LiCSBAS has limited utilization because it relies on sufficient LiCSAR data 
[19], [20]. 

However, some limitations have become essential matters. In this study, multi-temporal InSAR 
processing using LiCSBAS succeeded in detecting surface deformation at 100 landslide locations out of a total 
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of 358 landslide locations. The 100 landslide locations, on average, have an area of more than | hectare with 
an open landslide area without vegetation cover. The 100 landslide locations, on average, have an area of more 
than 1 hectare with an open landslide area without vegetation cover. This condition is related to the LiCSBAS 
output deformation pixels of 100 meters [19]. In addition, LiCSBAS utilizes LiCSAR products derived from 
Sentinel-1 with C-band [19], [20]. In densely vegetated areas, the Sentinel-1 C-band data produces 
interferograms with low coherence. LiCSBAS will automatically eliminate interferogram pixels with 
coherence values of less than 0.5. Some researchers have shown that C-band data is more suitable for measuring 
surface deformation in urban areas [19], [20], [50]. In addition, the hilly morphology in the research area can 
cause errors such as shadows. However, the combination of the LiCSBAS module and overlay of lithological 
formation and slopes in the Pacitan study area succeeded in identifying surface deformation at 100 landslide 
locations. 


4. CONCLUSION 

This study presents the use of LiCSBAS, an open-source analysis module integrated with LiCSAR, 
as a novel approach to address the issue of surface deformation in landslide locations with various physical and 
environmental conditions. They are not in an urban area or a large area; instead, they are in the Pacitan Regency. 
In addition, landslide distribution was spatially mapped based on lithological formations and slope classes and 
the effect of vegetation cover dynamics, using superimposed and Leaf Area Index analysis. 

Based on the results of the overlay analysis, 280 locations of landslides, or 73.5 percent, occurred in 
volcanic rocks and volcanic sediments, especially from the Arjosari, Mandalika, and Jaten formations, which 
form steep slopes and are prone to landslide — these rock types are spread throughout the district, except in the 
western part. Meanwhile, as many as 219 landslide locations, or 61 percent, are in the 15-30° slope class. The 
frequency of landslide locations increases until the slope class is 15-30°, then decreases with increasing slope 
value. 
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